*** Set working directory ***
cd "C:\..."

*********************** Load Data **********************************************
clear
graph drop _all
cap drop all

import excel "Data.xlsx", sheet("data") firstrow

************************** Setup ***********************************************
* Choose impulse response horizon
local hmax = 36

* Time variable
gen monthly_date = mofd(date)
format monthly_date %tm
drop date
rename monthly_date date
sort date
tsset date
egen time = group(date)
drop if date == .

* Standardize MP Shock series - unit standard deviation
egen mps_std = sd(chg_sveny02)
gen mps = chg_sveny02/mps_std

sort date
tsset date

* Normalize GSCPI for 1998-2019 sample
rename gscpi gscpi_unscl
egen gscpi_mean = mean(gscpi_unscl)
gen gscpi_dm = gscpi_unscl - gscpi_mean
egen gscpi_std = sd(gscpi_dm)
gen gscpi = gscpi_dm/gscpi_std

* GSCPI mp interaction terms
gen gscpi_mps = L.gscpi*mps
gen gscpi_mps1 = L.gscpi*L.mps
gen gscpi_mps2 = L.gscpi*L2.mps

* Put variables in log form
gen cpi_inf = ln(cpi)
gen rgdp_growth = ln(ip)
rename unrt unrt_chg
rename retail retail2
gen retail = ln(retail2)
rename pce pce2
gen pce = ln(pce2)
rename core core2
gen core = ln(core2)

sort date
tsset date

* LHS variables
forvalues h = 0/`hmax' {
	gen rgdp_growth_`h' = f`h'.rgdp_growth
	gen cpi_inf_`h' = f`h'.cpi_inf
	gen unrt_chg_`h' = f`h'.unrt_chg
	gen retail_`h' = f`h'.retail
	gen pce_`h' = f`h'.pce
	gen core_`h' = f`h'.core
}
************************** LP regressions **************************************
/* Run the LPs */
*ssc install lincomest

* PCE
eststo clear
cap drop b1 b2 u1a u2a d1a d2a u1 u2 d1 d2 Years1 Zero1 b3 u3a d3a u3 d3 Years2 Zero2
gen Years1 = _n-1 if _n<=`hmax'
gen Zero1 =  0    if _n<=`hmax'
gen b1=0
gen b2=0
gen u1a=0
gen u2a=0
gen d1a=0
gen d2a=0
gen u1=0
gen u2=0
gen d1=0
gen d2=0
gen stde1 = 0 

gen Years2 = _n-1 if _n<=`hmax'
gen Zero2 =  0    if _n<=`hmax'
gen b3=0
gen u3a=0
gen d3a=0
gen u3=0
gen d3=0

forv h = 0/`hmax' {
	newey pce_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(1/2).pce l(0/2).unrt_chg l(0/2).retail l(0/2).ebp, lag(`h')
replace b1 = _b[mps]                    if _n == `h'+1
replace u1a = _b[mps] + 1.645* _se[mps]  if _n == `h'+1
replace d1a = _b[mps] - 1.645* _se[mps]  if _n == `h'+1
replace u1 = _b[mps] + 1* _se[mps]  if _n == `h'+1
replace d1 = _b[mps] - 1* _se[mps]  if _n == `h'+1

replace b2 = _b[mps] + _b[gscpi_mps]                    if _n == `h'+1
matrix C1=(1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
matrix Vcomb1 =C1*e(V)*C1'         
mat SEcomb1 =cholesky(Vcomb1)      
scalar temp1 = SEcomb1[1,1]        
replace stde1 = temp1            

replace u2a = _b[mps] + _b[gscpi_mps] + 1.645* stde1  if _n == `h'+1
replace d2a = _b[mps] + _b[gscpi_mps] - 1.645* stde1  if _n == `h'+1
replace u2 = _b[mps] + _b[gscpi_mps] + 1* stde1  if _n == `h'+1
replace d2 = _b[mps] + _b[gscpi_mps] - 1* stde1  if _n == `h'+1

replace b3 = _b[gscpi_mps]                    if _n == `h'+1
replace u3a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d3a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u3 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d3 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* Core PCE
eststo clear
cap drop b4 b5 u4a u5a d4a d5a u4 u5 d4 d5 Years3 Zero3 b6 u6a d6a u6 d6 Years4 Zero4
gen Years3 = _n-1 if _n<=`hmax'
gen Zero3 =  0    if _n<=`hmax'
gen b4=0
gen b5=0
gen u4a=0
gen u5a=0
gen d4a=0
gen d5a=0
gen u4=0
gen u5=0
gen d4=0
gen d5=0
gen stde2 = 0 

gen Years4 = _n-1 if _n<=`hmax'
gen Zero4 =  0    if _n<=`hmax'
gen b6=0
gen u6a=0
gen d6a=0
gen u6=0
gen d6=0

forv h = 0/`hmax' {
	newey core_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(1/2).core l(0/2).unrt_chg l(0/2).retail l(0/2).ebp, lag(`h')
replace b4 = _b[mps]                    if _n == `h'+1
replace u4a = _b[mps] + 1.645* _se[mps]  if _n == `h'+1
replace d4a = _b[mps] - 1.645* _se[mps]  if _n == `h'+1
replace u4 = _b[mps] + 1* _se[mps]  if _n == `h'+1
replace d4 = _b[mps] - 1* _se[mps]  if _n == `h'+1

replace b5 = _b[mps] + _b[gscpi_mps]                    if _n == `h'+1
matrix C2=(1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
matrix Vcomb2 =C2*e(V)*C2'         
mat SEcomb2 =cholesky(Vcomb2)      
scalar temp2 = SEcomb2[1,1]        
replace stde2 = temp2            
replace u5a = _b[mps] + _b[gscpi_mps] + 1.645* stde2  if _n == `h'+1
replace d5a = _b[mps] + _b[gscpi_mps] - 1.645* stde2  if _n == `h'+1
replace u5 = _b[mps] + _b[gscpi_mps] + 1* stde2  if _n == `h'+1
replace d5 = _b[mps] + _b[gscpi_mps] - 1* stde2  if _n == `h'+1

replace b6 = _b[gscpi_mps]                    if _n == `h'+1
replace u6a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d6a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u6 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d6 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}

************************* Plot IRF's *******************************************
*** PCE ***
* Both states
twoway ///
(rarea u1a d1a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u1 d1  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b1 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ip_avg, replace

twoway ///
(rarea u2a d2a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u2 d2  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b2 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Tight Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ip_high, replace

gr combine fig_ip_avg fig_ip_high, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("PCE", size(small)) ///

gr rename fig_ip_both, replace


* Differential
twoway ///
(rarea u3a d3a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u3 d3  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b3 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("PCE: Interaction Coefficient", size(small)) ///

gr rename fig_ip_int, replace

*** Core PCE ***
* Both states
twoway ///
(rarea u4a d4a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u4 d4  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b4 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_cpi_avg, replace

twoway ///
(rarea u5a d5a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u5 d5  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b5 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Tight Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_cpi_high, replace

gr combine fig_cpi_avg fig_cpi_high, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("Core PCE", size(small)) ///

gr rename fig_cpi_both, replace


* Differential
twoway ///
(rarea u6a d6a  Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u6 d6  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b6 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("Core PCE: Interaction Coefficient", size(small)) ///

gr rename fig_cpi_int, replace

* Combine all 4 IRFs
gr combine fig_ip_both fig_cpi_both fig_ip_int fig_cpi_int, ///
rows(2) cols(2) graphregion(color(white)) plotregion(color(white)) ///

********************************************************************************